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FERROMAGNETIC  MITIGATION  OF  ELECTROMAGNETIC 

PULSE  EFFECTS 


I .  BACKGROUND  AND  MOTIVATION 

Electronic  equipment  in  regions  up  to  2000  miles  in  diameter 
can  be  made  inoperative  by  electromagnetic  pulses  ( EMP )  generated 
during  nuclear  detonation.  The  EMP-induced  transients  produce 
electromagnetic  fields  having  electric  field  magnitudes  up  to  50 
kilovolts  per  meter.  Such  fields  can  produce  permanent  damage  or 
malfunction,  particularly  to  solid  state  devices.  As  a  result, 
entire  communications  systems  can  be  wiped  out  by  a  single 
nuclear  burst.  Moreover,  EMP  can  produce  currents  as  large  as 
thousands  of  amperes  in  space  and  ground  based  electronics 
systems.  Such  currents,  if  not  attenuated  properly,  can  destroy 
any  semiconductor  device. 

There  are  a  number  of  techniques  which  can  be  employed  for 
EMP  mitigation.  These  include  the  use  of  extension  shields, 
semiconductors  (e.g.,  zener  diodes),  spark  gaps,  thermal  and 
electromechanical  devices  and  filters.  Unfortunately,  none  of 
the  presently  used  devices  are  adequate  for  systems  requiring 
fast  rise  times  and  simultaneous  large  power  handling  capability. 

Ferromagnetic  shields  may  provide  a  solution  to  the  EMP  prob¬ 
lem.  They  can  be  designed  to  reflect  part  of  the  energy  and 
absorb  the  rest.  Parameters  associated  with  these  mechanisms  are 
adjustable  within  very  broad  ranges  by  selecting  off-the-shelf 
materials  and  appropriate  device  configurations.  Moreover, 
magnetic  media  specifically  tailored  for  optimizing  EMP  mitiga¬ 
tion  should  provide  even  further  improvements.  These  improve¬ 
ments  can  be  accomplished  by  tailoring  parameters  such  as 
magnetic  permeability,  domain  switching  time,  magnetic  relaxa¬ 
tion,  saturation  magnetization  and  electrical  conductivity.  New 
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materials  like  metglass  (an  amorphous  ferromagnetic  material)  can 
be  employed  as  a  substitute  to  those  used  previously. 


Earlier  work  has  shown  [1]  that  one  can  obtain  a  42  db 
enhancement  of  the  EMP  shielding  factor  in  off-the-shelf  ferro¬ 
magnetic  conductors  over  that  in  nonf er romagnetic  conductors. 
Alternately,  the  use  of  ferromagnetic  shielding  can  reduce  the 
weight  by  an  order  of  magnitude  over  that  required  for  nonferro¬ 
magnetic  conductors  having  the  same  EMP  shielding  factor. 

Although  this  earlier  work  clearly  demonstrates  the  utility  of 
ferromagnetic  shields  devised  from  off-the-shelf  materials,  a 
realistic  understanding  of  the  behavior  of  such  shields  is  indeed 
lacking.  This  deficiency  is  principally  due  to  the  nonlinear 
nature  of  the  problem  which  led  earlier  workers  [1-5]  to  analytic 
solutions  which  oversimplify  the  magnetization  processes 
involved.  Such  oversimplifications  include  models  which  assume 
constant  or  a  limiting  value  of  the  magnetic  permeability.  In 
contrast,  realistic  descriptions  which  incorporate  actual  materi¬ 
al  characteristics  may  be  obtained  with  the  use  of  computer 
solutions.  The  present  report  describes  such  a  computer  solution. 
In  our  approach,  we  have  developed  a  theoretical  model  which 
allows  for  the  inclusion  of  magnetic  relaxation  as  well  as  other 
physical  parameters. 

II.  THEORETICAL  FORMALISM 
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We  then  have  from  Maxwell's  equations 
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where  yu  and  a  are  the  permeability  and  conductivity,  respectively, 
of  the  magnetic  film.  At  this  time  we  also  assume  that  fj  is  time 
independent.  Using  the  geometry  of  Fig.  2  we  get 
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Finally,  from  Eqs.  la  and  2a,  we  find 
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Taking  the  Laplace  Transform  of  Eq.  3,  we  get 
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whose  solution  is  of  the  form 
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To  solve  for  A  and  B,  we  introduce  the  boundary  conditions  for 
the  electric  and  magnetic  fields  at  y  =  0  and  y  -  d. 
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where  Zq  is  the  free  space  characteristic  impedance, 
simple  algebraic  "manipulations"  we  obtain 
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Substituting  Eq.  (5)  into  Eqs .  (7a)  and  7b)  we  get 
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This  relationship  is  very  much  analogous  to  a  relationship  found 
in  circuit  analysis.  We  define  f(s)  as  the  transfer  function  of 
the  magneto-conductive  medium.  We  define 
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We  now  wish  to  put  f(s)  into  a  form  which  can  be  transformed  back 
into  the  time  domain.  Substituting  the  above  definitions  into 
Eq.  8  we  obtain,  after  some  manipulations 
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where  the  summation  comes  from  the  expansion  of  the  factor 
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which  appears  in  the  denominator  of  an  earlier  stage.  Here  we 
have  let 


and 


•1  Z0, 


c2  -  ad  . 


Finally,  by  expanding  all  the  factors  containing  /s",  we  get 
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From  reference  [6],  Table  29.3.87,  we  have  for  the  Laplace 
transform 


n-1 


-k/sT 


.th 


_k 

e  Tt  . 

_  ft  (  *  ) 

2n  /  n+1  n  2/t 

'nt 


where  H  is  the  n  order  Herraite  Polynominal 
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the  transform  of  f(s),  we  have 
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where  and  M2  are  the  numerical  coefficients  in  Eqs.  (10a)  and 
(10b),  respectively. 

Also,  we  define 

f(t)  -  fx(t)  -  f2(t)  . 

Now,  we  have  to  convolute  the  transfer  function,  f(t  with  the 
incident  pulse  to  obtain  the  transmitted  pulse,  i.e.. 
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This  completes  the  derivation. 
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III.  RESULTS  AND  DISCUSSION 

The  convolution  integral  developed  in  the  last  section  ( Eq . 

12)  has  been  evaluated  for  three  input  pulses  (a),  a  rectangular 
pulse,  (b)  a  triangular  pulse,  and  (c),  a  double  exponential 
EMP  pulse  [7],  Each  of  these  cases  will  be  discussed  in  detail 
below  following  a  description  of  the  computer  program  used  to 
evaluate  the  convolution  integral. 

A  listing  of  the  computer  code  is  found  in  the  Appendix.  It 
consists  of  a  main  program,  CONVEMP,  and  two  function  subrou¬ 
tines,  FLAP  and  FLAP2.  FLAP  calculates  the  transfer  function 
f(t-r)  and  FLAP 2  calculates  the  pulse  h(t).  There  are  three 
versions  of  FLAP 2  listed,  one  for  each  of  the  three  pulses 
mentioned  above.  In  addition  to  these  user  subroutines,  several 
library  routines  are  called  to  evaluate  the  integral  and  plot  the 
results.  D01AHF  is  a  NAG  [8]  Library  subroutine  which  integrates 
the  function  FLAP 2  between  the  limits  TMIN  and  TB  to  the  relative 
accuracy  EPSR.  In  our  program,  FLAP 2  calls  on  FLAP  to  produce  the 
function  FLAP2  -  FLAP 2  *  FLAP,  and  it  is  the  product  function 
which  is  integrated.  The  integral  is  evaluated  repeatedly  from  t 

-  0  until  some  t  when  the  intensity  of  the  attenuated  pulse 

max 

approaches  zero.  EZXY  and  ANOTAT  are  NCAR  [ 9 ]  Library  routines 
which  plot  K  points  of  the  function  FP  vs.  T  and  are  used  in  our 
program  to  plot  the  results  of  the  above  integrations. 

A  comment  is  necessary  about  the  summations  which  appear  in 
Eqs .  10  and  12.  Our  program  does  not  contain  the  sura  over  1  and 
m  explicitly,  but  includes  all  terms  up  to  1  +  m  -  1.  (Note  that 
the  s-dependent  part  of  Eq.  10,  the  part  that  undergoes  the 
transform,  is  dependent  only  on  the  sum  1  +  m) .  Higher  terms  in 
the  expansion  are  dependent  on  higher  inverse  powers  of  t  and 
have  a  small  effect  near  the  origin  and  can  be  neglected.  The 
sum  over  n  is  included  explicitly  and  for  most  cases  there  was  no 
significant  improvement  in  increasing  n  from  10  to  100.  However, 
for  certain  values  of  the  parameters  and  at  a  sufficiently  large 
t  there  was  a  small  error  which  could  be  eliminated  by  increasing 
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n  to  40.  For  all  the  cases  discussed  below,  a  value  of  n  -  10 
was  used.  Finally,  the  high  precision  (REAL*16)  used  in  FLAP  was 
needed  because  many  of  the  terms  are  products  of  factors  contain¬ 
ing  very  large  exponentials  which  exceed  the  range  of  lower 
precision  variables. 

The  last  step  (before  discussing  the  three  pulses)  is  the 
evaluation  of  the  constants  c^  and  C2  in  terms  of  the  material 
parameters.  We  list  the  appropriate  definitions  and  values 
below: 


and  therefore 
and 

Furthermore , 

exponentials,  2nc2  +  ay  and  2nc2  +  2C2  -  ay,  become  equal  to 
(2n  +  l)c2.  This  makes  both  exponentials  the  same  and  greatly 
simplifies  the  calculation. 

A  plot  of  the  transfer  function  appears  in  Fig.  3.  The  curve 

4  -9 

reaches  a  maximum  of  1.48  x  10  at  t  -  3.2  x  10  sec.  and  falls 

—  8 

to  14.6  at  t  -  3.0  x  10  sec. 
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when  y  =  d,  the  two  factors  which  appear  in  the 
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Fig.  3  Plot  of  Transfer  Function  vs.  Time. 
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26x10“ 


Now  let  us  discuss  the  three  cases.  The  rectangular  and 

triangular  pulses  represent  simple  pulse  shapes  chosen  to  check 

out  our  program.  The  EMP  pulse  (see  Eq.  13  below)  is  one  usually 

used  in  the  literature  [7].  The  rise  time  is  related  to  the 

235  , 

"unfolding  time"  of  the  fission  of  U  and  the  decay  time  is 
related  to  the  duration  of  the  gamma  rays.  The  widespread  use  of 
this  curve  makes  it  a  suitable  choice  for  the  presentation  of  our 
calculations . 

a)  Rectangular  Pulse 

The  incident  pulse  has  the  value  h  =  1  for  0<t<  1.0  x  10  7 
and  h  -  0  for  t>1.0  x  10-7.  Figure  4  shows  the  attenuated  pulse 

plotted  out  to  t  -  2.0  x  10  .  The  pulse  rises  to  within  one 

r  —  5  —  8 

percent  of  its  maximum  value  of  9.0  x  10  at  t  -  2.0  x  10  and 

falls  to  8.0  x  10“7  at  t  -  1.2  x  10-7. 

b)  Triangular  Pulse 

The  parameters  for  the  incident  pulse  are  h  -  0  at  t  -  0, 
h  -  1.0  at  t  -  1.0  x  10~8  and  h  =  0  at  t  >  1.0  x  10-7.  Figure  5 
shows  the  attenuated  pulse  plotted  out  to  t  -  2.0  x  10  7. 

As  in  the  case  a)  above,  the  slope  of  the  attenuated  curve  lags 
behind  that  of  the  incident  curve.  The  pulse  rises  to  a  value  of 

8.25  x  10-8  at  t  »  2.0  x  10-8  and  falls  to  3.3  x  10-8  at 

t  -  1.2  x  10“7. 

c)  EMP  Pulse 

The  incident  pulse  has  the  form  [7] 

h  «  A(e"at-  e“8t)  (13) 

with  a  =  4.0  x  108  and  8  -  476.0  x  108  (t  in  seconds).  FLAP 2 
contains  a  routine  to  normalize  the  peak  to  unity.  The  maximum 
occurs  at  t  =>  1.0  x  10  .  Figure  6  shows  the  attenuated  pulse 

plotted  out  to  t  -  1.0  x  10"6.  The  pulse  rises  to 

8.8  x  10-8  at  t  -  2.5  x  10-8  and  falls  to  1.6  x  10  8  at 
t  -  1.0  x  10~8. 
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C0PLO1 


TIME  (SEC) 

Fig.  4  Plot  of  h^/lK  vs.  Time  for  a  Rectangular  Pulse  having 

an  amplitude  of  1  for  0  Z  T  Z  1.0  x  10“7  sec.  and  0  for 
T  >  1.0  x  IQ'7  sec. 
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C0PLO2 


TIME  (SEC) 

Fig.  5  Plot  of  h-p/h.  vs.  Time  for  a  Triangular  Pulse  having  an 

amplitude  of  1  for  T  s  1.0  x  10“®  sec.  and  0  at  T  =  0  and 

-7 

T  >  1.0  x  10  sec . 
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C0PLO5 


Fig.  6  Plot  of  h^/h.  vs.  Time  for  a  pulse  having  the  form 

=  A( e“aT  -  e‘ST)  with  a  =  4.0  x  106,  fi  =  476.0  x  106 

and  A  chosen  so  that  h„,„„  =  1. 

max 
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In  all  of  these  examples  the  transmitted  pulse  is  not  only 
attenuated  but  delayed  and  distorted.  This  is  most  easily  seen 
in  the  transfer  function  itself,  which  shows  the  transmission  of 
a  delta  function  input  pulse.  The  distortion  results  from  the 
progressively  greater  attenuation  of  higher  frequency  components 
in -the  signal.  The  large  attenuations . obtained  result  primarily 
from  reflection  at  the  front  surface,  so  that  even  very  thin 
material  can  be  an  effective  shield  against  radiation  fields. 
With  the  method  of  calculation  used,  any  influence  of  the  back 
surface  on  this  reflection  is  included. 


As  shown  here,  we  have  developed  a  realistic  description  of 
EMP  shielding  by  ferromagnetic  materials,  including  the 
attenuation  of  three  pulses  by  a  typical  ferromagnetic  shield. 
Moreover,  the  form  of  this  description  readily  allows  for  the 
inclusion  of  specific  properties  of  ferromagnetic  materials 
(e.g.,  magnetic  relaxation,  domain  structure  and  magneto¬ 
striction)  in  calculations  of  shielding  effectiveness. 
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PROGRAM  CONVEMP 
CHARACTER  IDENT*8,YN*1,DUM*1 

DIMENSION  FP(2001),T(2001),DT(3),XLAB(2),IW(102),TN(10),NPI(10) 
REAL*8  Y,Y0,FLAP2,D01AHF,EPSABS,EPSREL,ABSERR,W(800),TMIN,TI, 
1FL(10),EPSR,RELERR,TB 
REAL*16  PI,PIR2,XMl,F0 
EXTERNAL  FLAP2 

COMMON  /FBLK/C1,C2,PIR2^LPHA,F0,NMAX,/FBLK3/TMIN,TI,FL12, 
1NI,NPI,TN,FL,NF 

DATA  Cl/8.6637E-9/,C2/1.9090E-4/,XLAB/’TIME7r/,YLAB/’F17 
CALL  DATE(DT) 

TYPE*, ’Cl  =  ’,C1,’  C2  =  \C2 
5  TYPE*, ’RUN  IDENTIFICATION  (FILENAME)?’ 

IDENT=’ 

ACCEPT  136,IDENT 
C  EXPANSION  OF  DENOMINATOR 

TYPE*, ’HIGHEST  VALUE  OF  N  IN  SUMMATION?’ 

ACCEPT*, NMAX 

TYPE*,’  INPUT  FOR  FLAP2’ 

NF=0 

C  ACCEPT  INPUT  FROM  PULSE  SUBROUTINE 
Y0=FLAP2(0.) 

PRINT  100,CHAR(27),CHAR(15) 

PRINT  105, DT 
PRINT  145,IDENT 
PRINT  160, NI 
PRINT  115, NMAX 
NF=-1 

C  PRINT  INPUT  FROM  PULSE  SUBROUTINE 
Y0=FLAP2(0.) 

PRINT  130,CHAR(12),CHAR(18),CHAR(27) 

NP=0 

TYPE*, ’SUPPRESS  PRINTING  (Y/N)?’ 

ACCEPT  140, YN 
IF(YN.EQ.’Y’)  NP=1 
TMAX=TMAX*1.E— 6 
IDENT(7:7)=’*’ 

XM1=-1. 

PI=QACOS(XMl) 

PIR2=SQRT(PI) 

T(1)=0. 

IPT=IMAX+1 

Y=0. 

FP(1)=0. 

EPSR=l.E-5 

TMIN=0. 

NLIMIT=0 

IFAIL=0 

IT1=1 

T00=0. 

IIT=0 
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C  EVALUATE  CONVOLUTION  INTEGRAL  AT  POINT  TI  IN  INTEGRATION  REGION  I 
C  NEXT  I 

DO  20  1=1, NI 

T0=(TN(I+1)-TN(I))/NPI(I) 

IPT=NPI(I)+1 

IF(I.GT.l)  IIT=IIT+NPI(I-1) 

C  NEXT  TI 

DO  20  IT=2,IPT 
-  TI=TN(I)+T0*(IT— 1) 

K=IIT+IT 

T(K)=TI 

ALPHA=-C2 

Y=0. 

NCJ=0 

C  PERFORM  INTEGRATION  OVER  EACH  INTEGRATION  REGION  J  UP  TO  POINT  TI 
DO  30  J=1,NI 
NF=J 

IF(NCJ.EQ.l)  GO  TO  30 
IF(FL(J+1).LT.0.)  GO  TO  30 
TMIN=TN(J) 

TB=TN(J+1) 

IF(TI.LE.TN(J+1))  THEN 

TB=TI 

NCJ=1 

ENDIF 

TYPE*, TO, TN  (I),TI,TB 
C  NAG  INTEGRATION  ROUTINE 

Y0=  D01AHF(TMIN,TB,EPSR,NPTS,RELERR,FLAP2,NLIMIT,IFAIL) 

Y=Y+Y0 

30  CONTINUE 
20  FP(K)=Y 

C  NCAR  PLOTTING  ROUTINES 

CALL  ANOTAT(XLAB,YLAB,1,0,0,%REF(DUM» 

CALL  EZXY(T,FP,K,%REF(IDENT))  ' 

IF(NP.EQ.O)  THEN 

PRINT  100,CHAR(27),CHAR(15) 

PRINT  120 

C  PRINT  OUTPUT  AT  SELECTED  POINTS  ONLY 
DO  10  J=1,K,10 
IF(J.EQ.l)GO  TO  10 

IF(MOD((J— 1),40).EQ.0)  PRINT  125,(T(I),FP(I),I=(J-30),J,10) 

10  CONTINUE 

PRINT  130,CHAR(12),CHAR(18),CHAR(27) 

ENDIF 

TYPE*, ’ENTER  ANOTHER  SET  OF  DATA  (Y/N)?’ 

ACCEPT  140, YN 
IF(YN.EQ.’Y’)  GO  TO  5 
100  FORMAT(lX,Al,’I5i’,Al) 

105  FORMAT(’  DATE  =  ’2A4,A1) 

115  FORMAT!’  LAST  TERM  IN  EXPANSION,  N  =  ’,15) 

120  FORMAT(5X,4(4X,’TIME’,10X,’FCONV’,7X)) 
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125  F0RMAT(5X,4(E12.4,2X,E14.5,2X)) 

130  FORMAT(lX,3Al,’(4i’) 

135  FORMAT(A8) 

140  FORMAT(Al) 

145  FORMAT(lXrA8) 

160  FORMATC  NO.  OF  INTEGRATION  REGIONS  =  ’,15) 

STOP 

END 

C  THIS  ROUTINE  CALCULATES  THE  TRANSFER  FUNCTION 
REAL*8  FUNCTION  FLAP(TAU) 

REAL*16  CK1VX1,EX1,F,F1,F2,F3,F10,F1S,F2S,F3S,PIR2 
REAL *8  TI,FLAP2,TAU,TI0,TMIN 

COMMON  /FBLK/Cl,C2,PIR2rALPHA,F,NMAX,/FBLK2/Xl,EXl,CKl,FlS,F2S,F3S, 
1/FBLK3/TMIN,TI0 
TI=TI0— TAU 

IF((TI**2).EQ.O..OR.TI.LT.O.)  THEN 
FLAP=0. 

RETURN 

ENDIF 

F=0. 

F1S=0. 

F2S=0. 

F3S=0. 

CKl=l./(PIR2*TI**1.5) 

ALPHA=-C2 
DO  10  N=0,NMAX 
ALPHA=ALPHA+2.*C2 
ALP  H  A2=ALP  H  A*ALP  HA 
X1=ALPHA2/(4.*TI) 

EX1=EXP(-X1) 

F10=C1*CK1*EX1 

F1=-2.*F10 

F2=C2*F10*(C2+6*C1)*(2?N+1)**2/TI 
F3=-C1*C2**3*F10*(2*N+1)**4/TI**2 
F=F+F1+F2+F3 
F1S=F1S+F1 
F2S=F2S+F2 
F3S=F3S+F3 
10  CONTINUE 
FLAP=F 
RETURN 
END 


18 


C  THIS  ROUTINE  GENERATES  A  RECTANGULAR  PULSE 
REAL*8  FUNCTION  FLAP2(TAU) 

REAL *8  TAU,TB,T,FL(10),TMIN,FLAP 
DIMENSION  TN(10),NPI(10) 

COMMON  /FBLK3/TMIN,TB,FL12,NI,NPI,TN,FL,NF 
DATA  TN(l)/0./,FL(l)/l./ 

C  INPUT  PULSE  PARAMETERS 
IF(NF.EQ.O)  THEN 

TYPE*, ’NUMBER  OF  INTEGRATION  REGIONS?’ 

ACCEPT*, NI 

TYPE*, ’FOR  EACH  REGION,  NUMBER  OF  POINTS,  VALUE  OF  FUNCTION, 
1AND  END  POINT?’ 

DO  104=1, NI 

TYPE*, ’REGION  NO.  \I 

ACCEPT*, NPI(I),FL(I+1),TN(I+1) 

10  TN(I+l)=TN(I+l)*l.E-6 

RETURN 
ENDIF 

C  PRINT  PULSE  PARAMETERS 
IF(NF.EQ.-1)  THEN 
PRINT  110 

PRINT  115,(NPI(I),TN(I+1),FL(I+1),I=1,NI) 

RETURN 

ENDIF 

C  CALCULATE  PULSE 
FLAP2=FL(2) 

FL12=FLAP2 

FLAP2=FLAP2*FLAP  (TAU) 

110  FORMAT(2X,’NPTS,,4X,’TMAX’,9X,’FL’) 

115  F0RMAT(1X,I5,2E12.3) 

RETURN 

END 

C  THIS  ROUTINE  GENERATES  A  TRIANGULAR  PULSE 
REAL*8  FUNCTION  FLAP2(TAU) 

REAL *8  TAU, TB,T,FL(10),TMIN, FLAP 
DIMENSION  TN(10),NPI(10) 

COMMON  /FBLK3/TMIN,TB,FL12,NI,NPI,TN,FL,NF 
DATA  TN(l)/0./,FL(l)/l./ 

C  INPUT  PULSE  PARAMETERS 
IF(NF.EQ.O)  THEN 

TYPE*, ’NUMBER  OF  INTEGRATION  REGIONS?’ 

ACCEPT*, NI 

TYPE*, ’FOR  EACH  REGION-NUMBER  OF  POINTS, VALUE  OF  FUNCTION, 
1AND  END  POINT(MICROSECONDS)?’ 

DO  5,I=1,NI 

TYPE*, ’REGION  NO.  ’,1 
ACCEPT*,  NPI(I),FL(I+1),TN(I+1) 

5  TN(I+l)=TN(I+l)*l.E-6 

RETURN 
ENDIF 
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C  PRINT  PULSE  PARAMETERS 
IF(NF.EQ.-1)  THEN 
PRINT  110 

PRINT  115,(NPI(I),TN(I+1),FL(I+1),I=1,NI) 

RETURN 

ENDIF 

C  CALCULATE  PULSE  AT  POINT  TAU  IN  REGION  NF 
GO  TO  (10,20)  NF 
10  A1=FL(2)/TN(2) 

FLAP2=A1*TAU 

FL12=FLAP2 

FLAP2=FLAP2*FLAP(TAU) 

RETURN 

20  DELT=TN(3)-TN(2) 

A2=-FL(2)/DELT 

B2=FL(2)*TN(3)/DELT 

FLAP2=A2*TAU+B2 

FL12=FLAP2 

FLAP2=FLAP2*FLAP  (TAU) 

110  FORMAT(2X,’NPTS’,4X,’TMAX’,9X,’FL’) 

115  F0RMAT(1X,I5,2E12.3) 

RETURN 

END 

C  THIS  ROUTINE  GENERATES  A  DOUBLE  EXPONENTIAL  EMP  PULSE 
REAL*8  FUNCTION  FLAP2(TAU) 

REAL *8  TAU, TB,T,FL(10),TMIN, FLAP, GAMMA, GAMMA1, EX, C1,A 
DIMENSION  TN(10),NPI(10) 

COMMON  /FBLK3/TMIN,TB,FL12,NI,NPI1TN,FL,NF 
DATA  TN(l)/0./,TN(2)/l.E-8/,TN(3)/l.E— 7/,NI/3/,A0/l./ 

SAVE  A 

C  INPUT  PULSE  PARAMETERS 
IF(NF.EQ.0)  THEN 

TYPE*, ’NUMBER  OF  POINTS,  END  POINT  (MICROSECONDS)?’ 
ACCEPT*, NPI(4),TN  (4) 

TN(4)=TN(4)*l.E-6 
TYPE*, ’ALP  HA,  BETA?’ 

ACCEPT*,  FL(1),FL(2) 

DO  10,1=1,2 
10  FL(I)=FL(I)*1.E6 

NPI(1)=25*NPI(4)/100 

NPI(2)=NPI(1) 

NPI(3)=2*NPI(1) 

C  NORMALIZE  CURVE  TO  A0 
GAMMA=FL(2)/FL(1) 

GAMMA1 =GAMMA- 1 
TMAX=LOG(GAMMA)/(FL(l)*GAMMAl) 

EX=1.+1./GAMMA1 

C1=GAMMA1/GAMMA**EX 

A=A0/C1 

RETURN 

ENDIF 
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C  PRINT  PULSE  PARAMETERS 
IF(NF.EQ.-1)  THEN 
PRINT  105,NPI(4),TN(2) 

PRINT  100,FL(1),FL(2) 

PRINT  110,C1,TMAX 

RETURN 

ENDIF 

C  CALCULATE  PULSE  AT  POINT  TAU 

FLAP2=A*(EXP(-FL(1)*TAU)-EXP(-FL(2)*TAU)) 

FL12=FLAP2 

FLAP2=FLAP2*FLAP  (TA  U ) 

100  FORMAT!’  ALPHA  =  \E12.3,’  BETA  =  ’,E12.3) 

105  FORMAT!’  NPTS  =  ’,15,’  T  =  \E12.3) 

110  FORMAT!’  E(MAX)/A  =  ’,F7.5,’  TMAX  =  ’.E12.3) 
RETURN 
END 
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